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A novel sign-free Monte Carlo method for the Hubbard model has recently been proposed by 
Corney and Drummond. High precision measurements on small clusters show that ground state 
correlation functions are not correctly reproduced. We argue that the origin of this mismatch lies in 
the fact that the low temperature density matrix does not have the symmetries of the Hamiltonian. 
Here we show that supplementing the algorithm with symmetry projection schemes provides reliable 
and accurate estimates of ground state properties. 
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I. INTRODUCTION 

The understanding of low-temperature properties of 
doped Mott insulators is a central challenge in solid state 
physics. From the numerical point of view, this problem 
has remained elusive due to the sign problem. Recently 
Corney and Drummond have proposed a stochastic 
method in which the sign-problem does not explicitly oc- 
cur. They show that the density matrix of an arbitrary 
model Hamiltonian may be expressed as a positive sum 
over Gaussian operators. The imaginary time propaga- 
tion of the density matrix boils down to a Fokker-Planck 
equation governing the time evolution of the probabil- 
ity distribution in the space of Gaussian operators. One 
can then solve the Fokker-Planck equation by integrat- 
ing numerically the corresponding stochastic differential 
equation (SDE). For the Hubbard model on arbitrary 
lattice topologies and at arbitrary band fillings, the SDE 
has real stochastic and drift forces thereby leading to no 
explicit sign problem. 

The aim of this article is to test the precision of the 
method with respect to ground state properties. We will 
see - on the basis of simple examples - that the low tem- 
perature density matrix obtained by solving the SDE nu- 
merically does not have the symmetry of the Hamiltonian 
thereby producing biased ground state properties. The 
problem stems from the fact that a single Gaussian opera- 
tor breaks spin, lattice as well as translation symmetries. 
Since the weighted summation over the Gaussian opera- 
tors produces the density matrix, the summation has to 
restore the symmetries of the problem at hand. How- 
ever the summation is carried out stochastically, and it 
is a-priori not clear that the sampling is efficient enough 
to restore symmetries. At high temperatures this poses 
no further problems but as the temperature is lowered 
symmetry restoration fails. We show that one can solve 
this problem by projecting the density matrix onto the 
symmetry sector of the ground state. 

The organization of this paper is as follows. In the next 
section, we briefly review the formulation of the Gaussian 
Quantum Monte Carlo approach (GQMC) applied to the 



Hubbard model and in Appendix IU1 present the general- 
ization to four-fermion terms arising from Coulomb inter- 
actions. We will highlight the implicit assumptions used 
to derive the Fokker-Planck equation, absence of bound- 
ary terms, and discuss in Appendix IBI stochastic gauges 
which serve as a means to suppress them. In section ITTT1 
we show in detail how to implement the symmetry pro- 
jections and demonstrate their efficiency in section IIVI 
Finally we draw conclusions. 



II. THE GAUSSIAN QUANTUM MONTE 
CARLO METHOD FOR THE HUBBARD MODEL 

In this section we summarize the results of Ref. 0. 
Although the GQMC is general and can be general- 
ized to arbitrary Coulomb interactions as discussed in 
Appendix [21 we will concentrate here on the Hubbard 
model: 



H = c J Tc 
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St creates a fermion with 
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quantum numbers x = (i,cr) where i denotes the lattice 
site and a the z-componcnt of spin. Hence x runs from 
1 • • • N 8 = 2N, N being the number of lattice site. T is 
the hopping matrix. It is diagonal in spin indices and 
takes the value — t (— t') for next neighbors (next nearest 
neighbors). Finally a denotes a Pauli spin matrix and 

ct = (ct ct jj . Setting a — 1 yields the attractive Hub- 
bard model whereas setting a to a x or <j z the repulsive 
case [l2j . 

Corney and Drummond propose to expand the density 
matrix in terms of Gaussian operators: 



A(n) = det(l - n) : e 



(2) 



with n an N s x N s real matrix. The Gaussian operators 
are normalized, Tr A(n) = 1 and obey Wick's theorem 
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such that 

Tr[A(n)4c y = //., 

Tr A(n)cic y clc z 



Partial integration, under the assumption that boundary 
terms vanish, yields the Fokker-Planck equation for the 
probability distribution P(A,r): 



(1-nl .(3) 



The major result of Ref. 0] is that one can expand the 
density matrix in terms of a positive sum of Gaussian 
operators: 



(4) 



Clearly Tr [/5(r)] = Pi(t) grows exponentially with r. 
One can account for this exponential growth by attach- 
ing a weight factor to the Gaussian operators thereby 
obtaining: 



P(r) 



d\P{\, r)A(A), 



(5) 



with A = (fi,n), A (A) = f2A(n) and / dAP(A, r) = 1. 

The aim is now to formulate a stochastic process which 
samples the probability distribution P(A, r) in the space 
of Gaussian operators. To this end one considers the 
imaginary time evolution of the density matrix 



-P(r) 



1 r 
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so that in conjunction with Eq. (J5J we are left with the 
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The 



evaluation of the anti-commutator - 

anti-commutator can be transformed into a differential 
form acting on A(A): 
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The form of the diffusion matrices, D 



Pih,r) 
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(x,y),(w,z) 



A(A), 



T,r C x!y C w,z and D fx,y),{w,z) = Ei -^£,3/^^ 1S impor- 
tant. It depends on the manner in which we have written 
the Hubbard interaction term, or in other words on the 
choice of the diffusion gauge. The fact that the diffusion 
matrices factor out as above allows., us to fgrmulate the 
SDE. Furthermore, the fact that B l and C % are real for 
real values of n will lead to positiveness of the weights. 
The appropriate choice of diffusion gauges for general 
Hamiltonians is considered in Appendix IH1 

The assumption of vanishing boundary terms is essen- 
tial to justify the approach and boils down to the re- 
quirement that the probability P(A, r) has tails decaying 
sufficiently fast as |A| — * 00. At this point, one has to 
recall that the Gaussian basis is overcomplete, such that 
different probability distributions -P(A, r) will yield the 
the same density matrix. This degree of freedom is re- 
flected in a stochastic gauge invariance which is reviewed 
in Appendix [5] Hence, even if boundary terms appear 
the hope remains of eliminating them by an appropriate 
stochastic gauge choice. 

To proceed, let us assume that we can neglect the 
boundary terms. The Fokker-Planck equation can con- 
veniently be transformed into an Ito SDE 0], 

dfl = -flh(n)dT, (11) 
dn = -Adr + ^B (?) dVF ? + ^C (?) dW|, 



h{n) = Tr(A(n)H 



A = i n (T -UM)n + ^n(T -UM)n, 
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(8) with Wiener increments {dWf) = (dWl) = (dWjdW4) = 
0, and (dWldWt) = (dWjdW ? ) = dr% ? . Eq. (TTJ de- 
scribes the time evolution of walkers in the space of Gaus- 
sian operators. At r = 0, p(r = 0) oc 1 such that all the 
Walkers can be parameterized by A = (1, 1/2). At imag- 
inary time r they are distributed according to P(X, t) so 
that we have access to the density matrix. In particular, 

(9) any equal time observable is given by: 



In the above, n = 1 — n and 
M, 7 



(6) 



E.Tr 


A(A,)0 


E.Tr 


A(A,)" 



(12) 



77,77' 



(V 2 - n )(l n ),(lr,') (J <y,r l ^r l ',a'- 



where the sum runs over the set of walkers generated 
by the SDE. Since Wick's theorem applies for a single 
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Gaussian operator the numerator of the above equation 
may easily be calculated. 

As apparent from Eq. (fTTjl the weight of a Walker at 
imaginary time r reads 



2x2, U/t-4, <n>-l 



-/J-dT'fe(n(r')) 



(13) 



Since the "equal time Green functions", n, are real, h(n) 
is real and the weight remains positive. Hence the algo- 
rithm shows no explicit manifestation of the sign prob- 
lem. However, the weights grow exponentially with imag- 
inary time thus yielding an exponential increase in the 
variance. To circumvent this problem, we have adopted 
the reconfiguration scheme proposed in Ref. Q. In this 
approach the population of walkers is kept constant. 
Walkers with large weights are cloned and those with 
small weights suppressed in such a way that in the large 
population limit the density matrix remains invariant. 
Finally, after reconfiguration the weights of all walkers is 
equal to their average. 

We can now test the accuracy of the method on a 2 x 2 
Hubbard model (See Fig.^J. As apparent from Fig.QJa) 
at high temperatures the GQMC result for the energy 
(bullets in Fig. da)) compares well with the exact re- 
sult (solid line). However at low temperatures there is a 
systematic deviation. We have failed to account for this 
mismatch by i) enhancing the number of walkers ii) using 
different schemes for the integration of the SDE iii) vary- 
ing the imaginary time step iv) setting a to a x instead 
of a z in Eq. and finally v) using different stochas- 
tic gauges (see Appendices iBl and lUll . Adding additional 
noise terms by means of diffusion gauges only seemed to 
make things worse. Concerning point ii) we tried implicit 
and explicit Euler schemes and a higher order Milstein 
integrator £|. As the latter changes the order of the 
algorithm from 0(N 3 ) to 0(N 4 ) without reducing the 
systematic error, we prefer the Euler schemes. 

To acquire more insight into the origin of the mis- 
match, we compute the charge and spin susceptibilities: 
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FIG. 1: a) Energy as a function of inverse temperature 
as obtained from exact diagonalization (solid line), from the 
GQMC (bullets) and from the GQMC supplemented by the 
symmetry projection (squares). Here we have projected onto 
the total spin s — state and d-wave lattice symmetry, (b) 
The charge susceptibility and (c) the longitudinal and trans- 
verse spin susceptibilities. In the above, we use 60'000 walk- 
ers, an imaginary time step of Art = 0.0001, an explicit Euler 
scheme with adaptive time step, and a = a z in Eq. (0 



where N = J^jclcr and S a = J^jclcr a c?. Those quanti- 
ties are plotted in Fig. IHb),(c). As apparent, the charge 
susceptibility as well as xt follow rather precisely the 
exact result. On the other hand, xV diverges as 1/T 
thereby signaling that the low temperature density ma- 
trix has non-vanishing overlaps with s > spin sectors. 
To solve this problem so as to produce accurate ground 
state results we propose to implement symmetry projec- 
tion schemes. 



III. SYMMETRY PROJECTIONS 

Here we will assume that the low temperature density 
matrix has a large overlap with the ground state density 
matrix and a small admixture of excited states. If this 
assumption is correct, then projection onto the symme- 
try sector of the ground state will filter out the excited 
states and produce an accurate estimate of low temper- 
ature properties. Let us note that symmetry projection 
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schemes have been used successfully in the framework of 
the path-integral renormalization group approach where 
the ground state wave function is approximated by a sum 
of Slater determinants [f| . Here, we first review the math- 
ematics of symmetry projections and then show how to 
implement them in the context of the GQMC. 

Let us first consider finite groups with elements R and 
irreducible representations T> a (R). Group theory then 
tells us that 



(15) 



where l a corresponds to the dimension of the representa- 
tion. For continuous groups, the sum has to be replaced 
by the invariant integral: J2 R — > J dR 0. 

To show how symmetry projections rely on the above 
identity let us first consider the group of translations by 
lattice vectors R. 



f{R)£f{R)^=cl s 



i+R,a' 



(16) 



with f(R) = e <fl ' E *" pc f,^" and ct = -4= V-e ipt c- 

In the above, N denotes the number of lattice sites. 
Since the group of translations is an Abelian group, the 
irreducible representations are one-dimensional and la- 
beled by the total momentum K. Classifying states 
in Fock space according to their total momentum, K, 
yields T> 6 (R) = (i£,ag\f(R)\K,ag) = e liiit . Here 
1 = YIk a 1-^' a ii\' wnere a g labels all the 

states in Fock space with total momentum K. The pro- 
jection operator onto the Hubert space with total mo- 
mentum Kq reads: 



P K 



^Y,^o\f(R)\K o yf(R). 



This expression may readily be verified: 



(17) 



(18) 



6 K,K 



E jj^(Ko\f(R)\K )HK,a R \f(R)\K,a g ) 
<{K,ag\*)\K,ag) = ^2{K ,ag o \^)\K ,ag o ). 



Within the very same framework, we can define the 
projection on the Hubert space with total spin s. We first 
parameterize the rotations in terms of the Euler angles, 
w = (a, (3, 7), such that with 



a spinor transforms as: 



(19) 



(20) 



Here, S z corresponds to the total z-component of 
spin, J^i h. < '~ ,<jZc v an< ^ a s i muar definition holds for 
S y . Using Eq. (fTHf) and noting that V s mm ,(uj) = 

(s, m\T(u})\s, m!) where the quantum numbers m,m' de- 
note the z-component of spin, the projection onto the 
Hilbert space with definite spin s, and vanishing z- 
component of spin reads: 



2s + 1 
7d7 



du)(s,0\f{oj)\s,oyf(u)). (21) 



Since we have chosen to parameterize rotations in terms 
of Euler angles the invariant integral reads: f dui = 
J Q 2n da §q d/3sin(/3) J Q 27T d"f and 



(s,0|f(w)| S ,0> = P s (cos(/3)), 



(22) 



where P s denotes the s th Legendre polynomial. 

Since the GQMC method is a grand canonical ap- 
proach we have equally implemented projection onto 
fixed particle number. To this purpose, we define the 
gauge transformation: 



f{(j>) = e ^ S ' c h 



(23) 



such that T(4>)clT : (0) = e^ct. Projection onto a given 
particle number sector then reads: 



Pn = ^- [ (N\f(<f>)\N)if '(0). 
2tt Jo 



(24) 



Finally, we have implemented the C4 lattice symme- 
tries to classify states according to i) s-wave: even under 
parity and ir/2 rotations ii) d-wave: even under parity 
and odd under tt/2 rotations and hi) p x + ip y : odd un- 
der parity and acquires a phase factor e m l 2 under 7r/2 
rotations. We denote this projection by Piatt- 

Since the Hubbard Hamiltonian is invariant under lat- 
tice vector translations, spin rotations, gauge transfor- 
mations, 7r/2 rotations, the ground state density matrix 
will have definite momentum, spin, particle number and 
lattice symmetry. Our aim is now to project the density 
matrix produced by the GQMC onto a given symmetry 
sector and then use the projected density matrix, 



PpP\ 



(25) 



to compute observables. Here, P is a product of all or 
only some of the above symmetry projectors with general 
form, 



(26) 



P = / dxg(x)T(x), 



where T is unitary and P> = P. 

To simplify the calculation, we will assume that the 
observable O commutes with P: 



P,0 



0. 



(27) 
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2 x 2, 17/* = 4 
(n) = l 


GQMC + Symm. Proj. 
s — 0, d-wave 


Exact 


Energy/t 

5(7T,7r) 

iV(7r,7r) 


-2.1021 ±0.0007 
2.1933 ±0.0010 
0.2667 ± 0.0004 


-2.1026 
2.1947 
0.2664 



C//t = 4,t'/t = 0, (n) = 1 



TABLE I: GQMC with symmetry projection for the 2 x 2 
half-filled Hubbard model at U/t = 4. Here we have projected 
onto the d-wave and spin-singlet Hilbert spaces. To impose 
the spin projection we have to integrate over the three Euler 
angles. This integration is done numerically by replacing the 
three dimensional integral by a Riemann sum over 5 3 points. 
The thus produced systematic error is not included in the 
error bars. The results and error bars stem from averaging 
the data over imaginary time (squares in Fig. 0. 



such that 



(0>J 



Tr 


PpPO 


Tr 


PpO 


Tr 


PpP 


Tr 


Pp 



(28) 



since P 2 = P. Estimating the right hand side of the 
above equation boils down to the calculation of Pp where 
p si A(A) and the sum runs over the walkers pro- 
duced by integrating the SDE. Hence, using the result of 
Appendix El 



(O)p = 





f(x)A(X)6 


Y,xJdxg(x)Tr 


'f(x)A(X) 



J2 x fdxg(x)Tr A(A(x))6 



(29) 



where T(x)A(X) = A(A(as)). 

We test the above procedure on the 2x2 lattice of 
Fig. ^ As seen in Fig. HJa) (solid squares) by project- 
ing onto the spin-singlet and d-wave state we obtain a 
very accurate estimate of the ground state energy already 
at pt = 5. Averaging over subsequent imaginary times 
yields the results presented in Table [I] It is important to 
note that not only the ground state energy is very well 
reproduced but that also reliable estimates for the spin 
and charge structure factors, 



(30) 



N ® = ^£ e? ' (? ~ J) ' 



n->- n->), 



1,3 



are obtained. 



IV. ACCURACY TESTS 

Here we provide further tests triggered at assessing 
the accuracy of the method. We first consider systems 



L = 4 


GQMC + Symm. Proj. 
■3 — n P — n 


Exact 


Energy/t 
S(7r,7r) 

N(ir tt) 


-13.630 ±0.016 
3.66 ±0.013 
386 ± 001 


-13.6224 
3.64 
0.385 


L = 6 


GQMC ± Symm. Proj. 
s = 0, P = 


PQMC 


Energy/t 
S(7r,7r) 
iV(7r,7r) 


-30.87 ±0.04 

5.86 ±0.05 
0.400 ± 0.004 


-30.87 ±0.02 

5.82 ± 0.03 
0.418 ±0.025 



TABLE II: Comparison between GQMC and benchmark re- 
sults for the 4x4 and 6x6 Hubbard model. For both param- 
eter sets, we project onto total spin s = and total momen- 
tum P = 0. The L = 4 (L — 6) simulations were carried out 
with 12'000 (6'000) walkers, an explicit Euler scheme and an 
imaginary time step Art = 0.0005 (Art = 0.001). The exact 
diagonalization results for the L — 4 lattice stem from Ref. 
0- For the L = 6 lattice we compare with the auxiliary field 
projector QMC (PQMC) algorithm. 



where the sign problem is absent in auxiliary field QMC 
methods, that is, the particle-hole symmetric Hubbard 
model. Table |H] presents results at half-filling for both 
4x4 and 6x6 lattices. In both cases, one sees that the 
agreement with benchmark results (exact diagonalization 
for the 4x4 lattice and auxiliary field projector QMC 
(PQMC) for the 6x6 lattice) is excellent. Furthermore, 
the real space spin-spin correlations agree very well with 
the benchmark results (see Fig. [2J. 



U/t = i,t'/t = 
(n) = 0.625 


GQMC + Symm. Proj. 
s = 0, s-wave, iV = 10 


Exact 


Energy/t 
S(n,n) 
N(n,n) 


-19.576 ±0.012585 
0.737 ± 0.002 
0.5075 ± 0.001 


-19.584 
0.73 


U/t = 8,t'/t = -0.3 
(n)=l 


GQMC + Symm. Proj. 
s = 0, P = 0, s-wave 


Exact 


Energy/t 
S{n,n) 
N(n,n) 


-8.498 ±0.012 

5.09 ±0.07 
0.191 ±0.004 


-8.4884 
4.985 
0.1920 


U/t = 8,t'/t = -0.3 
(n) = 0.875 


GQMC + Symm. Proj. 
s = 0, P = 0, N = 14 


Exact 


Energy/t 
S(n,n) 
N(n,n) 


-12.01 ±0.40 
0.941 ±0.17 
0.266 ±0.01 


-12.50293 
0.964776 
0.27962 



TABLE III: Comparison between GQMC and exact diago- 
nalization results. Here we have used 12'000 walkers and a 
time step of Art — 0.0005. The GQMC is a grand-canonical 
simulation. Hence in cases where charge fluctuations are not 
negligible we project onto fixed particle number Hilbert spaces 
so as to allow comparison with exact diagonalization results. 
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4x4, U/t=4, <n>=l 



4x4, U/t=8, t'/t=-0.3, <n>=l 




(0,0) (1,0) (2,0) (2,1) (2,2) (1,1) 
r 

6x6, U/t=4, <n>=l 



Aux. Field PQMC 
GQMC + Syram. Proj. 





(0,0) (1,0) (2,0) (2,1) (2,2) (1, 

4x4, U/t=8, t'/t=-0.3, <n>=0.875 




(0,0) 



(0,3) 



(3,3) (0,0) 



'(3,0) (1,0) (2,0) (2,1) (2,2) (1,1) 



FIG. 2: Real space spin-spin correlations as obtained from 
the GQMC and comparison with benchmark results. See cap- 
tion of Table UTI for details of the GQMC simulations. 



FIG. 3: Real space spin-spin correlations as obtained from 
the GQMC and comparison with exact-diagonalization re- 
sults. See caption of Table ITTT1 for details of the GQMC sim- 
ulations. 



The crucial point is to show that in situations where 
the sign problem plagues the auxiliary field QMC, the 
Gaussian approach remains accurate. Table ITTll presents 
three data sets where the sign problem in the auxiliary 
field approach varies from mild to very severe. 

i) Let us start with the 4x4 Hubbard model with near- 
est neighbor hopping, t, and (n) = 10/16. The agreement 
between the GQMC and exact diagonalization is excel- 
lent. It is worth pointing out that in this specific case, the 
GQMC results with and without symmetry projections 
are identical meaning that the GQMC automatically pro- 
duces the ground state density matrix with correct sym- 
metries. We believe that this is due to the fact that at 
this large doping away from half-filling the ground state 
is very well described by a paramagnetic mean field so- 
lution. Such a mean field solution is exactly reproduced 
by the GQMC approach. 

ii) Our second example is the half filled 4x4 frustrated 
Hubbard model at U/t = 8. Here frustration stems from 
a next-nearest neighbor hopping t 1 jt — —0.3. Both Ta- 
ble lllll and Fig.[3fa) show that we obtain excellent agree- 
ment with the exact result. Note that for those model pa- 
rameters the finite temperature auxiliary field approach 
has an average sign of (sign) w 0.2 at (3t — 10 and of 



(sign) w 0.1 at (it = 15. 

iii) We now consider a parameter set which is out of reach 
for the auxiliary field approach, U/t = 8, (n) = 0.875 and 
t 1 jt = —0.3. Tablc llHl shows that we are capable of repro- 
ducing the exact results. However the fluctuations and 
hence the error bars in the MC data are large in com- 
parison to the half-filled case. Those large fluctuations 
stem from the symmetry projection. In particular, the 



Pf> 



is small and has large 



denominator in Eq. I|28|l , Tr 

relative fluctuations. In other words, the low tempera- 
ture density matrix (here we have propagated the walk- 
ers up to fit — 40) produced by the GQMC still includes 
many excited states, and it is hard to filter out ground 
state properties by imposing symmetries |13|. Neverthe- 
less, and as seen in Fig. I^b), comparisons with exact 
diagonalization results show that we are capable of ac- 
curately reproducing the details of real space spin-spin 
correlation function. 
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CONCLUSIONS 



tities of the Grassmann algebra |8|]: 



We have shown that the GQMC method produces in- 
accurate ground state properties since the numerical so- 
lution of the SDE fails to produce a low temperature den- 
sity matrix with the symmetry properties of the Hamil- 
tonian. To repair this sampling problem, we propose to 
a posteriori project the density matrix onto the sym- 
metry sector of the ground state. We have shown am- 
ple non-trivial tests, including situations where auxiliary 
field methods fail due to the sign problem, where this 
approach yields accurate and reliable results. Those re- 
sults confirm the point of view that the low temperature 
density matrix produced by the GQMC has a good over- 
lap with the exact zero temperature density matrix, but 
that the GQMC density matrix contains excited states, 
because the symmetries are not correctly reproduced. 
Those excited states are filtered out by the projection. 

There are many open questions which deserve further 
work. In particular, is it possible to improve the sampling 
by incorporating aspects of the symmetry projections di- 
rectly into the SDE? Also, we have not yet fully exploited 
the flexibility of the stochastic gauges. It is at present 
not clear if, with a suitable choice of stochastic gauge, 
the here mentioned symmetry problems may be solved. 
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APPENDIX A: UNITARY TRANSFORMATION 
OF A GAUSSIAN OPERATOR 

In this appendix we show that: 

e ld ' ht A(\) = A (A), (Al) 

with 



(£10 = e^W^', 
(£|:A(ct,c):|0 = A(?,Z')eM, 

i = /n d ^ e ^io<e(A2) 

J X 

Here £ x are Grassmann variables and |£) fermion coherent 
states. 

In a first step it is convenient to transform e ic hc into 
a normal ordered form. Since h is hermitian, h = UDU^ 
with D a diagonal and U unitary. With the canonical 
transformation 7^ = c'U we obtain: 

e *M = ]J ( » = JJ [1 + _ 1} ^t % ] 

X X 

= TT : e (° W '"- 1 )^" :=: ^ V'-ty. . 

= :e & (< " (A3) 

We can now compute the quantity e %c ^ hc : e ctBc : where 
B is an arbitrary matrix: 

gic+ftc . gC+Bc ._. e c+(e ih -l)c .. ^ Be ,__ 

1^1:^:17X71 = 

VZVr)Vie-? 6 -J*-^^\Z)e? &ihr 'e , >^ B+1 ^{7\ = 



mm >< 



(rj\ : e <=V(B+i)-i] C : | 7 )( 7 | =: e ^[e^ {B+i)-i]c . (A4) 

Here, we have carried out the substitution r) = e ih "q, bar- 
ing in mind that e lh is unitary matrix. 
The result of Eq. IA1I follows from: 



fidet(l-n) ..e-^i+e^iM^-ir^c.- 
det(l -n) A ^ t , _ „_ e t[ 2+(#i T_ ir i] g 



det(l — n) ^ 



dct(l — n) : e~ 



(A5) 



A(n) 



(n 7 -!)- 1 = \(e th - l)n T + ll (n T - l)" 1 



fl = Qdet 



(e lh - l)n T + 1 



Here, = h, A = (fl,n), and A = (fi,n). 

Before showing the above, let us fist recall some iden- 



APPENDIX B: DRIFT GAUGES 

Since the Gaussian operator basis is overcomplete, 
there are many probability distributions P(X, r) which 
will result in the same density matrix. This degree of 
freedom on P(A, r) is reflected in the choice of stochastic 
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gauges. Clearly, the aim is to find a gauge which will 
suppress boundary terms which could potentially show 
up in the partial integration step required to obtain the 
Fokker-Planck equation. Here, we introduce drift gauges 
and then propose some first ideas on how to choose the 
appropriate gauge. 

To formulate stochastic gauge invariance, it is useful 
to introduce the index [i : ■ ■ • such that for example 
A M=0 = n and A M = n x ^ for /z : 1 • • ■ N%. Then Eq. 
may conveniently be written as: 

-M^,A(A)1 

Z L J -\- 

+ Iyb®b®— 

~ o — M — v 




(Bl) 



- ^cfc® — 



A(A), 



with A = {Qh{n),A), B w = (0,B (l) ) and with C {l) = 

(0,C^). The above equation remains invariant under 
the transformation: 



g® = ( ,C (?) ) - 
A = (Qh(n),A) 



(Qg®,B & ), 
.(nfe(n),A + X; 



(B2) 



®B® +f®C (T> ), 



where g® and f® are arbitrary functions of n. This 
invariance stems from the fact that 



^A(A)=A(A). 
For a given stochastic gauge the Ito SDE reads: 



(B3) 



dn = -n 



h(n)dr - 5 (?) dW ? - f &AW { 



An = — 



A + J29 &B& +f & C & 



(B4) 



As apparent, one can modify the drift force from A 
to A at the expense of adding noise in the weights VL. 
Since our aim is to suppress the potential occurrence of 
boundary terms, one can follow the idea of choosing the 
gauge such that A prohibits the walkers, n, of drifting to 
infinity. In other words A x y should have the same sign as 
Tix,y Fulfilling this requirement for each pair of indices 
x, y leads to equalities. But since we only have have 



2N ( g % and f % for i : 1 • • • N ) degrees of freedom we 
can only fulfill the above condition on average. Defining 
a scalar product: 



(jlj A*J — ^ ^ Tlx,yA Xt ' 



we require that 



(n,A) 



> 0. 



(B5) 



(B6) 



APPENDIX C: DIFFUSION GAUGES - 
POSITIVE-P REPRESENTATION FOR THE 
GENERAL ELECTRONIC STRUCTURE 
PROBLEM 

One can also take advantage of gauge degrees of free- 
dom by adding terms to the Hamiltonian which cancel 
each other (or are identically zero) and hence do not af- 
fect the observables. If these terms are of fourth and 
second order in the fermionic operators, they add a con- 
tribution to the diffusion part of the SDE, which is com- 
pensated in the drift part, while the equation for the 
weight £1 remains unaffected. In this appendix we will 
show how such a diffusion gauge allows one to generalize 
the positive- P representation of the fermionic Hubbard 
model to the general electronic structure problem, in- 
cluding arbitrary hybridization, Coulomb and exchange 
terms. 

Since electronic structure calculations have important 
applications in quantum chemistry and the correspond- 
ing Hamiltonian 



+ ^ Yj V ^l4aCl al Cl a >C ja (CI) 



ijklaa' 

has even been dubbed the "theory of everything" [jj, 
a simulation approach without uncontrolled approxima- 
tions or systematic errors is highly desired. In Eq. IjClfl . 
Ci a denotes the destruction operator for an electron 
with spin a in the orbital i, the hopping amplitude 
and fi the chemical potential. The four-fcrmion terms 
cluC^rCia-iCja arise from Coulomb interactions. 

The many-body problem (|C1I) can be mapped to a 
system of stochastic differential equations with positive 
weights, by an appropriate choice of diffusion gauge, 
which generalizes the gauge proposed by Corney and 
Drummond for the Hubbard model Q, Il0| . Similar to 
section II, we choose a basis of Gaussian operators pa- 
rameterized by A, which can be represented in this case 
by block-diagonal matrices with N% elements for spin up 
and spin down, respectively (N Q is the number of or- 
bitals). The Fokker-Planck equation (|10|l for the proba- 
bility distribution P(X, r) reads 



eh 



P(X,t) = L[P(\,t)}. 



(C2) 
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For a suitably chosen diffusion gauge, the operator L in 
Eq. l|C2f) takes the form 



L 



d 1 d d 

12 j^ Aa + o!2 7^ B <*j^ B p 

^ d\ a 2 ^ d\ a d\p 

a a/3 r 



with B a and C a real coefficients (we present here the 
formulation which leads to the smallest number of noise 
terms). In this case, the Monte Carlo sampling can be 
done by integrating the Stratonovich SDE 

d\ a ( T ) = A a (\)dT+B a {X)dW(T)+C a (X)dW'(T), (C4) 

where (dW(T)dW'(r)) = and (dW(r) 2 ) = (dW'(T) 2 ) = 
dr. 

Gauge degrees of freedom can be used to modify the 
form of the operator L in (|C2(1 . In Ref . pj , the identity 



(C5) 



was used to map the Hubbard model to a system of real 
SDE. Here we will show how the identity 



(C6) 



can be used to obtain real SDE and positive weights f2 
for the more general Hamiltonian (jClfl . First, we note 
that the latter can be written as 



H 



X ^Ujh i ,j a + ^ Wi^kintjaflkla' , (C7) 



ijo 



ijklaa' 



with [fiija-, hkia- 1 ] — 0. In particular, for a = a' only terms 
with i ^ k,l and j / k.l appear. The tn a correspond to 
the chemical potential. We define 



'tijfjflijfj , 



(C8) 



and using Eq. I|C6|I express the ij kla cx'-contribution in 



Eq. I|C1[1 in the form 



denotes the sign of W-JIZ, 



— (n. 



hkla'Y 



\WS a A 

Lj £ Ll (8 ij h ij<J + 8 kl n kla ,). (C9) 



Each term H m (to = ija or ijklaa') gives a contribution 

im™ 1 "' to the drift term and the contributions Bm yp , 
Cm" to the diffusion term of the (stochastic) differential 
equations (|C4|) for Q and n xyp . No diffusion term appears 
in the equation of motion for f2 and we can write 



dn 

d7 



dr 



12 A ^ 



(CIO) 



J2 (Kr p + B%vz m + c^Cn) , (en) 



where the £ m , £' m are independent Gaussian random vari- 
ables with variance 1/dr. The hopping term l|C8|) yields 
only the drift-contributions 

■^■ijcr = tij<T n ij<n (C12) 
-^■ija P = ~2~ [ n xja (8iy ~ Tliya) + (8 x j ~ n xja) n iya\8 pa ■ 

(C13) 

With the gauge choice l|C9|l . the interaction terms yield 
a Fokker-Planck equation of the form of Eqs. I|C2|I and 
(IC3j) with drift terms 

A ?jkiaa' = -W$°,n ija n k i0< +WX,n i ivn k j<r>5 <T < T ', (C14) 



A'l'xyp 

ijklaa' 



^^J^xja^yi ^iya^j H~ {^xj ^xja^^v 



X(rkjer - SUkla' - Sij/2)S pcr 
+ [n x ^'((5fc y — Thkya') + {5x1 ~ ^xla'^kya'] 

x(n k ia' - sriija - S k i/2)5 p(T 'Y ( C15 ) 
and diffusion terms 



ijklaa' 



c 



ijklaa' 



^ ^ kla' 1/^ ^xjai^iy Tliya)5pa 

w ij<T f n x i(j'(dky — nkya')5 pa ' , (C16) 

Tlxja^Tliyafipa 



-S w ijo (8 x l - n x la')n k ycr'Spcr' . (C17) 



Note that the right hand side of Eq. IjClOp is —h m (n)fl, 
with h m (n) = Tr(A(n)if m ). Furthermore, since the 



l ija 



are real variables, which remain real during the integra- 
tion of Eq. lfCTT|) . it follows from Eqs. (|CT2)l . lfC14|) and 
(jClOl! that the weight Vl will always stay positive. 

For the actual implementation, it is simpler to use the 
Ito SDE, which is also numerically more stable. In this 
case, the drift terms have to be modified as 



.4 



Ito 



A 



Stratonovich 







dX 

d 



■0 



(C18) 



while the diffusion terms remain unaffected. The four- 
fermion term l|C9|l thus yields a contribution 



.4 



Tl x yp ,Ito 



VV kla' 



^ [^a?jV {^iy Tliya^) 



ijklaa' 2 

H~ i^xj T^xja ^^iya^^kla' $pa 
+ [nxla'{$ky — Kkya') 

a' J ' L ija$ pa' 



+ (8x1 - n xl<r ')n ky , 

\\R'xlcr(fiiy ^iya") i^xl ^xla^^iya^^kja 
-\- {jh%j(j{&fay — TT'/cycr) H~ i^xj nxja)nkya)nila]3aa'3pa^ 

(C19) 
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to the right hand side of (|C11|1 . 

We have tested the Ito-SDE for the two-site model 

H = H t + Hfj + H u + H c + H x + H s + Hi n + Hh 2 , 

(C20) 

where the individual terms are 

H t = ~t^(h 12a +h 21a ), (C21) 

H n = ^(nna + n 2 2<,), (C22) 
H u = \- ^2(uihu a hu- a + u 2 n22 a n 2 2~cr) 1 (C23) 

£7 

H c = v c ^2{h lla h 22c j + un CT n 2 2-ff), (C24) 

#x = ^^"^cr^l-o-, (C25) 

i7 s = w s 2^(ni 2o -ni2- CT + n 2 ian 2 i-a), (C26) 

(7 

#/u = hx ^(nn a ni2-a + n lla h 2 i-a), (C27) 

(7 

#/v 2 = ^-2 ^(«22o-ni2- C r + "-22CTW21-0-)- (C28) 
cr 

The (stochastic) differential equations IjClOfl and I|C11|) 

were integrated using an implicit Euler scheme and adap- 



tive time steps. Every n roC oiif steps, the family of ran- 
dom walkers was reconfigured according to the method 
of Ref . Q . While the high-temperature behavior is cor- 
rectly reproduced, similar and even more severe prob- 
lems as those observed in the simulations of the Hubbard 
model - most notably systematic errors in the energy at 
low temperature - also plague the simulation of the two- 
site model with Coulomb interactions. We have not yet 
checked whether symmetry projections can be success- 
fully applied in this case. 

For realistic electronic structure calculations we pro- 
pose to start with a Hartree Fock calculation of the 
electronic structure problem and to use both the occu- 
pied and unoccupied Hartree Fock orbitals in a subse- 
quent quantum Monte Carlo simulation. This is the same 
Hamiltonian and basis set used in full-CI (configuration 
interaction) or coupled cluster methods (CCM) in quan- 
tum chemistry, but the algorithm described here would 
enable to study a larger number of basis functions than 
in a full-CI calculation. Using the Hartree Fock density 
matrix instead of a multiple of the unit matrix as initial 
density matrix p(0) has the advantage that the initial en- 
ergy is already that of the Hartree-Fock solution, which 
will then be lowered further by the projection in imagi- 
nary time. Already a projection over a short imaginary 
time (3 (which can now not be interpreted as inverse tem- 
perature) will give an energy lower than the Hartree-Fock 
solution. 
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